;+
; NAME:
;       FFORE_LBD
;
; PURPOSE:
;       Calculate the Galactic 8-um foreground fraction at (l,b,d).
;
; CATEGORY:
;       IRDC distance modeling
;
; CALLING SEQUENCE:
;       ffore = FFORE_LBD(l, b, d,[,/NONNATIVE][,/HYPERFINE]
;                         [,/FROMFIT][,MT=mt]) 
;
; INPUTS:
;       L -- Galactic Longitude [deg]
;       B -- Galactic Latitude  [deg]
;       D -- Heliocentric distance [pc] (Either scalar or vector)
;
; OPTIONAL INPUTS:
;       MT -- Model type (see BGPS_FFORE.pro) [Default: !IRDC.mt]
;
; KEYWORD PARAMETERS:
;       NONNATIVE -- Calculate ffore for a single distance (native
;                    usage is to calculate ffore as a function of d.
;       FROMFIT   -- For use with one of the specialized fitting
;                    routines -- ignore.
;
; OUTPUTS:
;       FFORE -- Foreground fraction at the distance(s) requested.
;
; OPTIONAL OUTPUTS:
;       NONE
;
; COMMON BLOCKS:
;       FFORE_BLOCK  -- Common block associated with the IRDC modeling code.
;
; MODIFICATION HISTORY:
;
;       Created:  08/20/12, TPEB -- Cut code from FUNCTION
;                                   IRDC_MORPH_FFORE (subroutine in
;                                   IRDC_MORPH.pro).
;
;-

FUNCTION FFORE_LBD, glon, glat, d, NONNATIVE=nonnative, $
                    HYPERFINE=hyperfine, FROMFIT=fromfit, MT=mt
  
  COMMON FFORE_BLOCK,n,rb3,R0,dist,R,Z,tau,f_data,sig_f,corr,farlist,do_Tfit,$
     lpstr,rho_hi,rho_h2,rho_star
  
  FORWARD_FUNCTION BGPS_FFORE_MODEL, BGPS_FFORE_MODEL2, BGPS_FFORE_SIMPLE, $
     BGPS_FFORE_SIMPLE2, BGPS_FFORE_LOOPIE, STAR_DENSITY
  RESOLVE_ROUTINE, 'BGPS_FFORE', /COMPILE_FULL_FILE, /EITHER, /NO_RECOMPILE
  
  ;; Get galactic params
  defsysv, '!MW', exists = exists
  IF ~ exists THEN galactic_params 
  defsysv, '!IRDC', exists = exists
  IF ~ exists THEN irdc_params, SUFF=suff
  
  fromfit = KEYWORD_SET( fromfit )
  IF n_elements(first) EQ 0 THEN first=0b
  
  IF KEYWORD_SET( nonnative ) THEN BEGIN 
     cdist = d
     d = dindgen(!MW.NBINS * 50.)*(!MW.BINSIZE / 10.) + (!MW.BINSTART / 10.)
  ENDIF
  dist = d
  
  ;; Here is where we specify the values for the model we want to use
  A  = !IRDC.A
  B  = !IRDC.B
  C  = !IRDC.C
  NG = !IRDC.NG
  GW = !IRDC.GW
  Td = !IRDC.Td
  mt = !IRDC.mt
  
  IF ~fromfit || first THEN BEGIN
     first = 0b
     ;; And before we can calculate FFORE, we must get (R,Z) and rho's
     omni_lbd2rz, glon, glat, d, r, z
     IF mt EQ 0 || mt EQ 2 THEN BEGIN
        rho_hi   = hi_density(R,Z)
        rho_h2   = h2_density(R,Z)
        rho_star = star_density(R,Z) 
     ENDIF
  ENDIF
  
  CASE mt OF
     0: fj = BGPS_FFORE_MODEL(1.d3,A,B,C,NG,rho)
     1: fj = BGPS_FFORE_SIMPLE(1.d3,A,B,C,rho)
     2: fj = BGPS_FFORE_MODEL2(1.d3,A,B,C,rho)
     3: fj = BGPS_FFORE_SIMPLE2(1.d3,A,B,C,NG,GW,rho)
     ELSE: message,"You've got a problem with model type.  RUM-TIME!!!"
  ENDCASE
  rho_total = TOTAL(rho)
  
  IF KEYWORD_SET( nonnative ) THEN BEGIN
     
     ilt = WHERE( d LE cdist, nilt)
     ffore = (nilt NE 0) ? A + (TOTAL(rho[ilt]) / rho_total) : A
     d = cdist
     RETURN, (ffore < 1.d)
     
  ENDIF ELSE BEGIN
     ffore = dblarr(n_elements(d))
     
     ;; Need to loop through distances because of integrating to a distance
     FOR i=1L, n_elements(d) DO $
        ffore[i-1] = A + (TOTAL(rho[lindgen(i)]) / rho_total)
     
     RETURN, (ffore < 1.d)
  ENDELSE
END
